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Abstract 

Direct verification of the existence of an infinite set of multicritical non-perturbative 
FPs (Fixed Points) for a single scalar field in two dimensions, is in practice 
well outside the capabilities of the present standard approximate non-perturbative 
methods. We apply a derivative expansion of the exact RG (Renormalization 
Group) equations in a form which allows the corresponding FP equations to appear 
as non-linear eigenvalue equations for the anomalous scaling dimension 77. At zeroth 
order, only continuum limits based on critical sine-Gordon models, are accessible. 
At second order in derivatives, we perform a general search over all rj > .02, 
finding the expected first ten FPs, and only these. For each of these we verify 
the correct relevant qualitative behaviour, and compute critical exponents, and the 
dimensions of up to the first ten lowest dimension operators. Depending on the 
quantity, our lowest order approximate description agrees with CFT (Conformal 
Field Theory) with an accuracy between 0.2% and 33%; this requires however that 
certain irrelevant operators that are total derivatives in the CFT are associated 
with ones that are not total derivatives in the scalar field theory. 
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Circumstantial evidence strongly suggests that there exists an infinite set of 
multicritical non-perturbative FPs for a single scalar field in two dimensions, 
corresponding to the universality classes of multicritical Ising models, equivalently 
to the diagonal invariants of the unitary minimal {p,p+l) conformal models with 
p = 3, 4, ■ ■ ■ [1] [2], however direct verification of these facts is in practice well outside 
the capabilities of the standard approximate non-perturbative methods: lattice 
Monte Carlo, resummations of weak or strong coupling perturbation theory and the 
epsilon expansion. (The impracticableness of the epsilon expansion for higher p is 
covered in ref. [3], implying similar difficulties in weak coupling perturbation theory, 
while lattice methods suffer from difficulties of locating and accurately computing 
the multicritical points in the at least p — 2 dimensional bare coupling constant 
space). In this letter we demonstrate directly that these multicritical points do 
exist - by an approximation scheme which is reviewed below. Our main motivation 
for the present letter is to show that the approximation scheme, which we presented 
and applied to three and four dimensional scalar field theory in refs. [5] [6] , is powerful 
enough, in fact to automatically uncover, and reliably describe, the expected much 
richer set of non-perturbative theories in two dimensions. Many of our operator 
dimensions are computed to at least the same accuracy as that of the most relevant 
operator - very much in contrast to all the above approximation methods, suggesting 
that our equations may also provide a reliable description well away from FPs. 
These facts are certainly encouraging for our ambition of developing a reliable and 
accurate analytic approximation method of general applicability to non-perturbative 
quantum field theory [4] -[6]. 

As mentioned in the abstract we use the approximation scheme carried to the lowest 
sensible order, to verify that each of the first ten multicritical points (p = 3, • • • , 12) 
have the correct qualitative behaviour - e.g. shape of potential, number of relevant 
directions, expected parities under (p —(p etc. - and compute critical exponents, 
and scaling dimensions of up to the first ten lowest dimension (integrated) operators. 
All of the ~ 100 quantities agree with CFT to (sometimes much better than) 33%, 
with a weak improvement in accuracy with increasing multicriticality. However, 
we find that some integrated irrelevant operators appear in the scalar field theory 
spectrum, which can only be identified with operators which are total derivatives 
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in the CFT. It appears to be an interesting challenge to understand why these 
integrated operators do not vanish. 

We start with a review of the method[5]. We wish to formulate quantum field theory 
in a way in which it is obvious that our approximations are renormalizable[4]. We 
do so by never needing to determine bare quantities, such as a bare action Saq [v^] , 
or equivalently counterterms. In this letter we discuss only quantum field theories 
with no low energy mass scale. (The generalisation to massive ones is in principle 
straightforward). If indeed such theories are independent of the cutoff Aq then the 
(renormalized) theories have scale invariance. As is well known this is achieved only 
by assigning scaling dimensions to the operators, which are different from their naive 
dimensions in general. Thus for a single component scalar field (fi{x) in D (euclidean) 
dimensions we have that its scaling dimension is = ^{D — 2 + rj), where r] is the 
anomalous scaling dimension. For example, this implies by dimensional analysis 
that propagators, as a function of momentum q, appear as 

<(p (pXxl/q'^-'^ . (1) 

To write the condition that the theories be scale invariant, and to formulate an 
efficacious approximation scheme, it is helpful to introduce a (low energy) scale^ A. 
Now by dimensions, (1) appears as 

<VV>^ X^/(?VA') (2) 

for some function /. We will require that A is introduced in such a way that / can 
be Taylor expanded for small Q'^/A^, because our approximation will follow from a 
momentum expansion (equivalently derivative expansion) by dropping terms beyond 
some maximum order. This requirement implies that the infrared singularity in (1) 
has been smoothly regulated. Therefore A is equivalent to a smooth infrared cutoff. 
We introduce it into the partition function Z[J] as C{q, A), satisfying C{q, A) — > 
as — > 0, by writing SAoi'p] ^ ^Aoi^p] + ^(f.C~^.(p, where only C depends on A. 
From this it is straightforward to write down a differential equation for Z[J] with 
respect to A. It turns out however to be helpful to transform this to a differential 

In the Wilson sense this is the intermediate cutoff at which one defines the effective 
action [7]. This connection is developed in ref.[4]. 
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equation for the Legendre effective action. (One important reason is because the 
integrals involved will converge, with our choice of C, only if the full self-energy is 
used - expanded as a power series to the prescribed maximum order). In terms of 
this, after expressing y?, C and q as dimensionless quantities^ using A, we find[5] 



(— + d^A^ + Aa-D)r[(^] 



POO 

Jo 



Ciq^) dq 



-1 



(3) 



(q, -q) 



Here t = ln(///A) with fi some arbitrary reference scale. The angle brackets refer to 
an average over all directions of the vector q. C is a normalization factor, introduced 
for convenience, by a numerical rescaling. A,^ — f-j^ counts the number of fields 
in a given vertex, and Ag counts the number of derivatives in a given vertex[5]. 

The requirement of scale invariance is now simply given by ^r[(y?] =0. Substituting 
this and a derivative expansion 

V[^] = J d^'x {V{cp, t) + lid^ipfKi^, t) + {d^^YH^{if, t) + {UiffH^i^, *) + •••} , 

(4) 

into (3) and expanding the RHS (Right Hand Side) to some maximum order in 
derivatives yields n coupled second order non-linear ordinary differential equations 
for the fixed point coefficient functions V{(p), K{(p), Hi{(p), H2{<~p), ■ ■ ., where n is 
the number of undetermined coefficient functions. As such, we expect at first sight 
there to be a 2n parameter set of solutions. In fact, generally there are only 
a discrete set of possible solutions. Of course only a discrete set of solutions is 
generally expected on physical grounds: they correspond to the possible massless 
continuum limits (continuous phase transitions) with the prescribed field content. 
Mathematically this arises because nearly all choices of BCs (Boundary Conditions) 
for the differential equations lead to solutions with singular behaviour at some finite 
real value of the field (p [5] [6]. Another way of seeing why only a discrete set of 
solutions is allowed is as follows. We take for simplicity a theory with Z2 symmetry 
(p <-> — and assume d^p ^ 0. (See later for what happens when these conditions are 
relaxed). In this case there are n BCs given by symmetry - namely that the first 



C must be chosen to scale correctly [5]. 
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derivative of the coefficient functions {V'{(p), K'{(p), etc.) vanish at <^ = 0, while 
a further n BCs are given by the behaviour of the coefficient functions for large (p 
which are determined, up to a proportionality constant, by dimensional analysis - 
assuming that in this limit A can be neglected: 

K{ip)(X(p^^-^^/'^^-^, ■■■ as ip^oo. (5) 

Since these 2n conditions'^ are imposed on a 2n parameter set, we again expect 
only a discrete set of solutions. Actually this is not the whole story because the 
parameter rj in (3) (with ^r[(^] = 0) must still be determined. Generally this is 
only possible if there is a reparametrization invariance of the RG equations. This 
turns the equations into non- linear eigenvalue equations for rj. (To see this, use the 
invariance to fix an extra condition) . Clearly it is important that such an invariance 
is preserved by the approximation scheme. We can do this if we choose C(q'^) to 
be homogeneous in q, i.e. C{q^) oc q'^'^, since then a reparametrization invariance of 
the equations exists according to the following (non-physical) scaling dimensions: 

[a^] = [g^] = l, [^] = K + D/2, 

6 

hence [V]=D, [K] = -2{k + 1), [Hi] = -D - 4{k + 1) , 

and this is clearly not affected by neglecting higher derivative terms. The value of 
K may be determined uniquely, by considerations of convergence, to be the smallest 
integer larger than D/2 — 1 [5]. 

The operator spectrum may be determined by linearization of (3) about the FP 
solutions. By separation of variables the perturbations are of the form SV{(p,t) = 
e{iJ./A)^ v{ifi), 5K{ip,t) = £(/x/A)^ /c((^), 5Hi{ipA) = e{p/A)^ hi{ip), etc, for some 
functions f, k, hi etc, where e is infinitessimal and A appears as an eigenvalue 
for the linearized equations. In fact A is quantized since we have again 2n + 1 
constraints: n from symmetry {v'{0) = 0, k'{0) = 0, • • • for even eigenfunctions, 
or v{0) = 0, k{0) = 0, • • • for odd eigenfunctions), n from dimensional analysis 
{v{(p) oc ip^^-^)/'^v ^ k{(f) oc (^(^~2-A)/d^-2^ • • • as — > OO, again providing A can be 
ignored in this limit), and one from a normalization constraint allowed by linearity. 
Sjjt,^ plays the role of an infinitessimal coupling constant of dimension A, conjugate 

^ If Z2 symmetry is dropped, (5) for ip — > ±00 provides all 2n conditions. 
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to an integrated operator of the form ^[(p\ = J d^x {v{(p) + ^k{(f){d^(ff' + •••}. Note 
that this imphes that the operator itself has dimension D — X. 

In interpreting the operator spectrum it is important to discard any redundant 
operators[8]. These are operators of the form[5] 

^[<p] = ld''xFMSn<P]/S<p{^) , (7) 

corresponding to reparametrizations of the effective action. They have no physical 
significance, and no well-defined scaling dimension since it depends on the details of 
the RG used[8]. Since we already have 2n + 1 constraints, the requirement (7) will 
overconstrain the problem leading to no solutions unless the redundant operators 
exist for special reasons (viz. symmetries). We know of two such operators. One 
is the operator corresponding to the symmetry (6), which thus has (with some 
arbitrary normalisation) Fx[v] = —x^dfj_(f{yi) — {k + D/2)<f>{x.), v{ip) = DV{(f) — 
{k + D/2)(fV'{(p), k{(f) = -2{k + l)K{(p) + D/2)(pK'{(p), ■ ■ •, even parity, and 
eigenvalue A = 0. And the other [8] corresponds to the (p translation symmetry of 
the unsealed RG equations. It has F^[(p\ = 1, v{(p) = V'{(p), k{(p) = K'{(p), ■ ■ •, and 
odd parity. By operating with Jd^x 5/5(p{yi) on (3), it can be seen that, to any 
order of the derivative expansion, X = d^, as expected on dimensional grounds. 

From now on we set D = 2. This implies k=1. To 0(9^) we drop all derivatives 
higher than second on the RHS of (3). Since the higher FP coefficient functions 
Hi, ■ ■ ■, then satisfy linear equations given by the LHS of (3) with ^r[(/?] = 0, they 
must vanish, otherwise they have singularities at = (for general values of ry) which 
is unacceptable. Thus the inverse operator in (3) is the same as that computed in 
ref.[5]. Performing the average, the q integral, choosing — 1/4, and matching both 
sides of (3) we find for V{(fi, t) and K{ip, t): 

' \v'-2V=-(l-^)pKiV") (8a) 



9t 2^ V 4 

and ^ + l^K' + rjK= (l - |) lp^{V") 



V"K" - AV"'K 
2{V"^-AK) 



(86) 



2KV"V"'^ + 13V" K'^ - AOKK'V" „ _ KV'V'"'^ - AKK'V" + VK'"^ ^ 



^ 2{V"'^ -AKY ^ {V"'^ - AKf 

12KK" + lli^'^ 2QK'V"V"' - UKV'"^ - UK'^ KV'"^ - K'V'V" + K'"^ 

24K{V"^ - 4K) ^ 6{V"^ - 4Kf (F"^ - 4Kf 
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and tan ^ is taken in the range < tan ^ < tt. Eqns. (8) hold true only if 

K>0 and V" > -2Vk , (10) 

for otherwise the integrals diverge at unphysical poles. These conditions are 
sufficient to ensure that the obvious physical stability requirements are satisfied. 

Consider now the 0(5°) case, where the RHS of (86) is also dropped. Then r] = 0, 
since otherwise the FP K solution is singular: K{(p) (xl/(p^. With = 0, K{(p) is not 
determined. We will here simply assume that K{(f) = 1 is a good approximation. 
The remaining FP equation is now 2V{(p) = Pi{V"). Since -Pi (a) is a positive 
monotonically decreasing function which diverges at a = —2, the 'potential' U (x) = 
— J^dy Pi^{2y) is bounded below with a single stationary point. If we think of V 
as position x, and (p as the time r, then we see that, apart from the Gaussian FP 
{V{(p) = tt/S), V{(f) has only periodic solutions, as follows from Newton's equation 
cPx/dr'^ = —dU/dx. These correspond to the semi-infinite line of critical sine- 
Gordon models (that is circular bosons with a radius tunable to any value larger 
than some minimum radius). It would take us too far from our present purpose 
however to flesh this out. Since = 0, eqn. (5) does not apply; in fact the unsealed 
V{ip) is nowhere independent of A and vanishes in the limit A ^ 0. (Incidentally, the 
dipipV'{ip) term in (8a) corresponds, in the above Newtonian analogy, to a negative 
friction term proportional to velocity and to the time. Physically it is then easy 
to understand that this term is responsible for generic singular behaviour and that 
the careful balancing act (5) is necessary if this is to be avoided.) 

Returning to the full 0(5^) equations, we substitute a differentiated (8a) into the 
V'" terms in (86) to turn (8) manifestly into a pair of second order differential 
equations. We analyse their FPs by the numerical method of relaxation [9]. We can 
find the approximate p — 3 (Ising model) solution and its eigenvalue r] = r]i by using 
a full set of BCs and supplying an initial guess (e.g. the D — 3 Ising solution[5]), 

6 



0.5 
0.0 
-0.5 

V'(0) 

-1.0 
-1.5 
-2.0 

0.01 0.10 1.00 

Fig.l. V'{0) versus rj. Note that no solution was found for rj > 1.08. The p = 3 • • • 12 
FPs correspond to the 10 places where V'{0) crosses the axis. 

however finding the other solutions this way is much harder. (More details of the 
numerics are given in the appendix). Anyway, we want to demonstrate that for 
^ 0, the eqns.(8) have no other FPs except the {p,p + 1) points. This is the 
expected result, by Zamolodchikov's c-theorem[10], if the FP can be the result of 
flow from a Landau Ginsburg Lagrangian based about the c = 1 ultraviolet Gaussian 
FP. In fact we will here strictly only show that no other connected FPs exist in the 
following sense: Starting from our approximate Ising FP, we relax the BCs by 
dropping the requirement V'{0) = and instead determine a unique solution as a 
function of rj by using relaxation while gradually changing rj away from rji. The 
result is plotted in fig.l for all r] > .02. Each time V'{0) crosses the axis we store 
the solution and use this as a guess for 'polishing' by relaxation with the full set of 
BCs. (A more conclusive demonstration that there are no other FPs would involve 
relaxing e.g. the K'{0) condition also, and searching a two dimensional parameter 
space). Similarly, using relaxed BCs and starting from the exactly known operators 
given earlier or below, we can determine all^ the even and odd parity eigenoperators 
and eigenvalues at any of the FPs. Here one must be careful that the normalisation 

^ Here we can be confident that all solutions are connected since we do not expect or 
find at any point a two dimensional linear space of solutions to the relaxed BCs. 
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condition can be satisfied at any value A. We chose v'(0)^ + v(0)^ = 1, dropping the 
requirement v'{0) = (v(0) = 0) for the even(odd) parity solutions. 

Some results are displayed in the table below. The dimensions and parities of the 
two most relevant operators are not displayed since the unit operator v{(f) = 1, 
k = hi = ■ ■ ■ = 0, has dimension and even parity for any FP and to any order in 
the derivative expansion, while the field itself v{(f) = k = hi = ■ ■ ■ = 0, always 
has dimension dip (= r]/2 in this case) and odd parity. (These follow easily from the 
vanishing of the RHS of (3) linearized about the FP). The correlation length critical 
exponent is given by = 1/ A where A is the largest even parity eigenvalue, excluding 
the unit operator [7]. All numbers were determined to an accuracy greater than 
that shown.^ Also shown in the table are the corresponding results expected from 
CFT[1][2][11]. These follow from the dimensions A„,^ = ^^^^2^^"+!)^'"^ 
parities Pn,m = (-)(^+i)(p+i)+(™+i)p of the scalar primary fields ^n,m{x), $i,i = 1 
and $2,2 = <f, where m, n are integers in the range 1 < n < p — 1, 1 < m < p. 
In addition there are an infinite number of Virasoro descendents, however all but 
L_iL_i^n,m = n^»n,m havc too high a dimension to correspond to the 0(5^) results 
in the table. Thus 2p — 3 relevant operators are expected {^k,ki k = 1, • • ■ ,p — 1 
and ^k+i,ki k = 1, • • • ,p — 2), of which p — 1 are even and p — 2 are odd. We 
checked that this is true for all our solutions. In fact, as one can see from the table, 
these CFT values and our results match unambiguously. As further evidence for 
the correspondence we note that if the bare Lagrangian is taken to be of Landau 
Ginsburg type with a potential which is an even polynomial in </? of degree 2p — 2 
[1] [2], then we might expect the FP Wilson potential to have p — 1 minima; It does. 
Examples are shown in fig. 2 and ref. [12]. Similarly, the v{<f) component of the 
operators turn out to have the maximum number of nodes expected (c.f. fig. 4). 

We do not find all the irrelevant operators expected from CFT, which is why there 
are blank entries in the table. For each p we obtain a cluster of solutions and then a 
gap in dimensions, much larger than expected from CFT. We report the p = 3, 4, 5 
cases in detail. For p = 3, we find the first operator, other than those in the table, is 
an even operator with A = 20 ± 2. (The estimated error is due to truncation in the 

^ Thus the small dimensions at the bottom of the table correspond to A's computed 
to better than 5 significant figures, and agreeing with CFT to 3 significant figures. 
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Fig. 3. Fixed point kinetic factors K{ip) for p = 3, 4, 5. As in fig. 2, the full solution 
is shown as a full (or for p = 5 dashed) line, while the dotted parts of the curve are 
given by the asymptotic expansion (using Ay and Ak)- Eqn.(6) has been used to 
scale the solutions to conventional normalisation K{Q) — 1, from the K"gy{l) = 
normalisation. (See appendix). 
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Fig.4. The v{lp) component of the 3^'^ - 6*^ eigenoperators at the tricritical Ising 
model fixed point {p — 4). As in figs.l and 2, the solutions have been scaled to 
K{0) = 1 normalisation. Furthermore the even(odd) operators have been normalised 
so that f (0) = 1 (^'(0) = 1). Note that the l'^^ operator thus has v{ip) = 1, and the 
2"'^ has v{ip) = (p. The 3'^'^ - 5^^ are identified at the bare level with the powers 
- in refs. [1] [2] (plus lower dimension counterterms of the same Z2 symmetry), 
while the 6*^*^ is one of the "puzzling operators" . 

asymptotic BCs). In the p = 4 spectrum we find only those operators in the table, in 
particular no sign of L_2L-2$2,2 (dimension and parity = 4^ ). In p = 5, we 
find an ll*'^ and a 12*^^ operator with = 2.61+ and 2.80~ respectively, and then 
nothing else out to at least A = 8, and in particular the primary $41 (A^ = 6"*") 
is missing. However, high dimension operators can correspond to operators with 
many derivatives, and therefore may not appear at our level of approximation. 

More puzzling are the integrated operators that are there, but (apparently) should 
not be because they only match total derivative operators in CFT (given as boldface 
entries in the table). Thus at p = 3, the 4*^ operator's A^ matches well with n^2,2i 
and not with any other operator - in particular the first odd operator that is not a 
total derivative has A = 6|. Similarly all the other irrelevant dimensions mentioned, 
match □ on the relevant CFT operators, with the 11**^(12*'^) operators in p = 5 
matching the values 2|^(3^ ). Only for the 7^^ operator in p = 4, and the ll*'^ 
and 12*^ in p = 5, would alternative assignments be possible within the accuracy 
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suggested by the table: to $3,1 (3+), $3,1 (2| ), and $4^2 (3| ) respectively. 
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Table 1. 77, z/, and the A^'s of the 10 lowest dimension operators (ordered by 
increasing dimension), for the first 10 multicritical points. The 0{d^) answer is 
shown in the first row and the associated exact CFT result in the second row. 
Worst determined number: rj for p = 4 (33% off). Best determined number: u for 
p = 12 (0.2% off). ?7's accuracy gradually improves from p = 4 to p = 12. 1/ is worst 
determined at p = 3 (13%) after which all are determined to error less than 2% 
and decreasing with increasing p. The worst determined operator dimension is the 
3""^ at p = 5 (25%) after which errors decrease with increasing p and/or increasing 
dimension. The dimensions of the first half of the relevant operators are always 
overestimated while those of the second half are always underestimated. 

These puzzling solutions could be discarded if they were redundant, but, as argued 
below (7), this requirement overconstrains. Indeed, assuming that the most general 
(Poincare invariant) redundant operator (7) is of the form F^[ip] = ax'^d^ip{x) + 
f (v?(x)) + • • •, where a is a constant, / some differentiable function, and the rest 
results in terms only of 0(9^) or greater in 7[</?], then it is easy to show that (for any 
dimension D) ^[ip] is a redundant operator only if v{ipi)V{ip2) — V{(pi)v{(p2) for all 
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pairs of stationary points <^i, <^2 of V- These relations are not even approximately 
satisfied by our puzzling operators. Finally, using the cases p = 3, 4, 5, the 'best fit', 
before the large gap in dimensions, can be summarised as follows: All the relevant 
operators appear, together with a set of irrelevant states whose A^'s match □ 
acting on all these operators but two: the first (n$i^i = 0) and the last (□$p_i^p_2 
is missing). Note that the 0{d^) spectrum has an unbroken oscillating pattern of 

Z2 parities (H 1 — ■ — I — ), which would have been broken for A > 2 had it also 

(or only) uncovered the irrelevant primaries here. 




Fig.5. The k{(p) component of the 3^*^ - 6*^^ eigenoperators at the tricritical Ising 
model fixed point (7? = 4). Normalisations are those of fig.4. Recall that the 1^* and 
2'^'^ operator have k{(p) = 0. 



We briefly mention directions for further research. Evidently, taking the derivative 
expansion to higher order would be informative. We should, for example, begin to 
see approximate degeneracies in the irrelevant part of the spectrum. Expanding the 
flow equations to second order in the perturbations would allow computation of the 
operator product (fusion) coefficients, while non-scalar couplings would allow access 
to operators with non-zero spin and e.g. an investigation of the approximately 
satisfied Virasoro symmetry. The full partial differential equations (8), analyzed 
e.g. again by relaxation, can describe for example flow (crossover) from one FP 
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to another, using the present analysis to provide the BCs; in particular one can 
investigate this way the form of the bare Landau Ginsburg Lagrangian. Finally it 
would be interesting to consider, by relaxation, 'flow' of the FPs with (general real) 
dimension D. All of the above suggestions would help in resolving the reason for 
our "puzzling operators" . 

Note added in proof. 

A natural interpretation of our puzzling operators would be that they are examples 
of "shadow operators" [13]. (I thank T. HoUowood for bringing this work to my 
attention). These are operators that exist in the continuum theory but completely 
decouple in the conformal invariant limit. If so, we should expect to find that, in 
our approximation, their fusion coefficients with CFT operators are very small. 
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Appendix A. Numerics. 

We sketch some of the less obvious points involved in the numerical solution of (8) 
at, and linearized about, the FPs. Note first the reasons for relaxation: it is not 
necessary to cast the eqns (8) in the form V"{(p) = • • • and K"{(p) = • • •; two point 
BCs and eigenvalue problems are easily incorporated; and it is efficient in the sort 
of searches illustrated in fig.l, since the previous solution is a good guess for the 
next [9]. It can also cope with stiffness [9], which is there because small perturbations 
from a solution result in singularities at some finite (p. Conditions (10) can be 
violated in the iterands. If this happened it worked best to back track and shorten 
the jumps between iterations. It is necessary to develop the asymptotic BCs (5) at 
least to next-to-leading order, introducing explicitly proportionality constants, say 
Ay and A^-Th.is is because we set these BCs at some finite point (p = pasy^ and 
imposing just (5) would imply, by (8), singular values for the second derivatives. 
In fact we developed the asymptotics yet one further order so as to check that the 
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truncation error was well under control. (pASY must lie within a window bounded 
below by points where the asymptotic series is not a sufficiently good approximation 
and above by points where roundoff error (due to increasing stiffness) prevents the 
relaxation program from converging. This window rapidly shrinks with increasing p 
because both effects (we believe) depend on the size of ~ ip^^"^. We found empirically 
that the upper bound is just above the last turning point in K. This turning point 
always exists, so we ensure we are automatically inside this window by setting 
fASY = 1 and normalising with K'J^^y{l) = 0, where Kasy{p>) is the asymptotic 
expansion of K. We checked that our results are completely insensitive to moderate 
changes in ^pasy ~ 1- We factor out the rapid increase in V in the asymptotic 
region (5), by writing V{ip) = [V{0) + Ay f'^'^]V (ip) , where the integer r ^ 2/?7, the 
multiplying factor being differentiated analytically. Numerical errors are then under 
control even for p — 12 where asymptotically V ~ p'^^^. For the eigenvalues A, it 
was enough to substitute the FP solutions into the linearized equations, using the 
same mesh, imposing the leading asymptotic behaviour on v, k at (p = Pasy In the 
positive [negative] parity search, we recorded those cases where v{l) changes sign 
but v'{0) [v{0)] does not, since sometimes these are due to passing over eigenvalues, 
however cases where both v'{0) [v{0)] and v{l) changed sign were due to k{0) [k'{0)] 
passing through zero. The most rapid asymptotic behaviour was factored out of 
both V and k. For the search we chose (pasy = 1, but for 'polishing' it was often 
necessary to choose (paay < 1 for convergence. We required as much convergence as 
possible before choosing (paay just inside the remaining wild oscillations, checking 
insensitivity to moderate reductions in this value. Finally next to leading order v 
and k asymptotics were used as a check particularly that no solutions were spurious. 
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